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ABSTRACT 

The silhouette cast by the horizon of the supermassive black hole in M87 can now be resolved with 
the emerging millimeter very-long baseline interferometry (VLBI) capability. Despite being '--^ 2 x 10^ 
times farther away than Sgr A* (the supermassive black hole at the center of the Milky- Way and 
the primary target for horizon-scale imaging), M87's much larger black hole mass results in a horizon 
angular scale roughly half that of Sgr A*'s, providing another practical target for direct imaging. 
However, unlike Sgr A*, M87 exhibits a powerful radio jet, providing an opportunity to study jet 
formation physics on horizon scales. We employ a simple, qualitatively correct force-free jet model to 
explore the expected high-resolution images of M87 at wavelengths of 1.3 mm and 0.87 mm (230 GHz 
and 345 GHz), for a variety of jet parameters. We show that future VLBI data will be able to constrain 
the size of the jet footprint, the jet coUimation rate, and the black hole spin. Polarization will further 
probe the structure of the jet's magnetic field and its effect on the emitting gas. Horizon-scale imaging 
of M87 and Sgr A* will enable for the first time the empirical exploration of the relationship between 
the mass and spin of a black hole and the characteristics of the gas inflow/outflow around it. 

- techniques: interferometric — submillimeter — galaxies: active 



Subject headings: black hole physics - 
— galaxies: jets 

1. INTRODUCTION 

With the advent of millimeter VLBI, for the first time 
it has become possible to resolve the horizon of a black 
hole (Doeleman et al. 2008a; Doeleman 2008). Most 
of the theoretical and observational efforts on this front 
have been directed towards imaging Sagittarius A* (Sgr 
A*), the radio source associated with the supermassive 
black hole at the center of the Milky Way (see, e.g., Fal- 
cke et al. 2000; Broderick & Loeb 2005, 2006a,b; Miyoshi 
et al. 2007; Doeleman et al. 2008a). This is because Sgr 
A* subtends the largest angle of any known black hole 
candidate (roughly 55 /xas in diameter). Nevertheless, 
there is at least one other promising target for millime- 
ter VLBI imaging, namely M87's supermassive black hole 
(hereafter labeled 'M87' for brevity). 

At a distance of 16Mpc and with a black hole mass 
of 3.4 X 10^ Mq, the apparent diameter of M87's horizon 
is approximately 22 /ias, appearing about half as large 
as Sgr A*. Like Sgr A*, M87 is many orders of mag- 
nitude under-luminous in comparison to its Eddington 
luminosity. However, M87 is different from Sgr A* in 
five important respects. 

First, M87 is located in the Northern sky, and thus 
is more amenable to imaging by telescope arrays com- 
prised of existing millimeter and submillimeter observa- 
tories, which exist primarily in the Northern hemisphere. 
Thus, observations of M87 will generally have lower at- 
mospheric columns and longer observing opportunities. 
Arrays of existing telescopes (see §2 for details), spanning 
projected baselines of 10'*km, are already sufficient to 
resolve the silhouette at 1.3 mm (230 GHz) and 0.87 mm 
(345 GHz) with beam sizes of 17/ias and 11 /las, respec- 
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tively, along the minor axis. The inclusion of new facili- 
ties over the next few years will improve baseline coverage 
substantially, resulting in a roughly circular beam with 
these resolutions. 

Second, M87's mass implies a dynamical timescale 
of roughly 5 hr, and a typical orbital timescale at the 
innermost-stable orbit of 2-18 days, depending upon the 
black hole spin. Thus, unlike Sgr A*, the structure of 
M87 may be treated as fixed during an entire day, one 
of the underlying assumptions of Earth aperture synthe- 
sis. This means that it will be possible to produce se- 
quences of images of dynamical events in M87, occuring 
over many days. 

Third, M87 exhibits a powerful jet, providing a unique 
opportunity to study the formation of relativistic jets on 
horizon scales. Realistic ah initio computations of jet for- 
mation are not presently possible, requiring jet simula- 
tions to make a variety of assumptions about the under- 
lying plasma processes. Examples include assumptions 
about the dissipation scale, cooling rate, jet mass-loading 
mechanism, acceleration of the observed non-thermal 
electrons and even the applicability of the magnetohy- 
drodynamic (MHD) prescription (De Villiers et al. 2003; 
McKinney & Gammie 2004; Komissarov 2005; McKin- 
ney 2006; Hawley & Krolik 2006; Tchekhovskoy et al. 
2008; Igumenshchev 2008; McKinney & Blandford 2008; 
Nishikawa et al. 2005). Thus, horizon-resolving images of 
M87's jet will critically inform our understanding of the 
formation, structure and dynamics of ultra-relativistic 
outflows. 

Fourth, unlike Sgr A*, one is primarily driven towards 
submillimeter wavelengths by resolution requirements. 
There is no evidence for an analogous interstellar scat- 
tering screen that would blur high resolution images of 
M87 (cf. Bower et al. 2006). Nor is there evidence that 
M87 is optically thick at or about millimeter wavelengths 
(cf. Broderick & Loeb 2006a, b). Thus, modeling and in- 
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terpreting images of M87 is likely to be somewhat more 
straightforward . 

Fifth, despite being under-luminous, M87 is similar to 
other radio-loud quasars. The same cannot be said of Sgr 
A*, which is nearly dormant and is visible only due to 
its extraordinary proximity. As such, the faint low-mass 
black hole in Sgr A* is much more representative of the 
population of black holes that are accessible to the Laser 
Interferometer Space Antenna (LISA '^). Careful compar- 
ison between horizon-scale imaging of Sgr A* and M87 
will be crucial as a touchstone for comparing the astro- 
physical properties of these two populations, and how 
the black hole mass and spin affect the electromagnetic 
and kinetic properties of the surrounding plasma. 

It is not surprising, therefore, that already high- 
resolution images of M87 has been aggressively pursued. 
The best known of these are the 7 mm VLBI images ob- 
tained by Junor et al. (1999) (and more recently by Ly 
et al. 2004; Walker et al. 2008). 3mm and 2mm images 
have also been reported (Krichbaum et al. 2006; Kovalev 
et al. 2007) . These represent the highest resolution obser- 
vations of a relativistic jet ever produced, probing down 
to roughly 5 Schwarzschild radii. Nevertheless, they are 
presently insufficient to directly address many outstand- 
ing questions in jet-formation theory. Among these ques- 
tions are the role of the black hole spin in jet production, 
the structure of the jet forming region, the importance 
of magnetic fields and the particle content of the jet it- 
self. However, all of these uncertainties may be settled by 
pushing to submillimeter wavelengths. In this sense, mil- 
limeter imaging promises to provide a qualitative change, 
not simply an incremental improvement. The reason for 
this is simple: the minimum relevant physical scale for 
a black hole of mass M is GM/c^ , and thus resolving 
this scale provides full access to the physical processes 
responsible for jet formation and accretion. In addition, 
the observation of the characteristic silhouette associated 
with the black hole settle the present ambiguity in in- 
dentifying its position relative to the observed radio core 
(Marscher et al. 2008). 

Motivated by the above arguments, we present in this 
paper predicted theoretical images of M87 for a class of 
force-free jet models which fit all existing observations. 
In particular, we assess the ability of possible millimeter 
& submillimeter VLBI experiments to distinguish differ- 
ent critical jet parameters. In §2 we review how a mil- 
limeter array can be created out of existing telescopes. 
In §3 we describe the jet-disk model we employ for M87. 
In §4 and §5 we present images and polarization maps, 
respectively, for a variety of jet models. Finally, §6 sum- 
marizes our main conclusions. 

2. POTENTIAL MILLIMETER ARRAYS 

A detailed discussion of existing telescopes that may 
serve as stations in a millimeter/submillimeter Very- 
Long Baseline Array (VLBA) can be found in Doeleman 
et al. (2008b). These facilities are: CARMA in Cedar 
Flat, California; the SMTO on Mount Graham, Arizona; 
the LMT on Sierra Negra, Mexico; APEX & ASTE, soon 
to be followed by ALMA, on the Atacama plane in Chile; 
the IRAM dish on Pico Veleta (PV), Spain; the IRAM 
Plateau de Burc (PdB) array in France; and a number of 

^ http://lisa.nasa.gov/ 
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Fig. 1. — Tracks in the u-v plane (representing the Fourier- 
transform coordinates of the sky) made as the Earth rotates by 
existing and forthcoming telescopes as seen by M87. Blue tracks 
are associated solely with the existing North American telescopes 
(SMA, JCMT, CARMA & SMTO), green includes the two Eu- 
ropean observatories (PdB & PV), orange includes the LMT and 
yellow includes telescopes in Chile (APEX & ALMA). The tracks 
are labelled according to H (CSO, SMA, JCMT), C (CARMA), 
S (SMTO), B (PdB), V (PV), L (LMT) and A (APEX, ASTE, 
ALMA). 

telescopes atop Mauna Kea, including the JCMT, SMA 
and CSO, that may be used independently or phased 
together to provide single Hawaiian station. The mil- 
limeter flux from M87 is a factor of 2-3 lower than Sgr 
A*, and thus telescope sensitivity will be a key factor 
in resolving small-scale structure. This may be achieved 
in a variety of ways, including phasing together the in- 
dividual telescopes in the CARMA, ALMA, SMA and 
PdB into single apertures at their respective locations. 
In addition, a shift to large-bandwidths will likely be re- 
quired. 

As mentioned in §1, M87 has a number of observa- 
tional advantages over Sgr A*, the primary one being 
its location in the Northern sky. This is evident in the 
baseline coverage provided by the previously mentioned 
telescopes, shown in Fig. 1. M87 is simultaneously vis- 
ible at the European stations (PV and PdB) and the 
North American stations (Hawaii, CARMA and SMTO). 
With the inclusion of the LMT, the baseline coverage is 
nearly complete within a narrow, roughly East- West el- 
lipse, providing a resolution in this direction of approxi- 
mately 17/xas and 11 /zas at 1.3mm and 0.87mm, respec- 
tively. Introducing a station in Chile results in a similar 
resolution along the North-South direction. 

Critical for Earth aperture synthesis, M87's substan- 
tially larger mass implies that it may be regarded as a 
stationary source, even on horizon scales, for timescales 
of order a day. The sparseness of the baseline coverage 
of Sgr A*, which varies over 15 min timescales, is due 
not only to the limited number of existing stations, but 
also the fact that visibilities taken throughout the night 
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Fig. 2. — Spectra of our canonical jet-disk model. The disk & 
jet components are shown by the short-dashed, red line and the 
long-dashed, blue line, respectively. The total spectrum is shown 
by the thick green line. Data points are collected from Birctta et al. 
(1991); Despringre et al. (1996); Perlman et al. (2001, 2007); Tan 
et al. (2008). The errors on the data points represent the instru- 
mental uncertainty, not the variability, and thus are an underesti- 
mate for the multi-epoch SED. The spectra of the other jet-disk 
models are nearly indistinguishable from the one shown. 

cannot generally be assumed to arise from the same in- 
trinsic brightness distribution. Thus, in the absence of 
additional telescopes, efforts to interpret features in the 
image of Sgr A* (such as a black hole silhouette, Mills 
1962), will need to be performed in Fourier space. In con- 
trast, M87's comparably long dynamical timescale means 
that images can be produced and analyzed directly. This 
is even true during dynamical periods in which the ac- 
cretion flow is changing as rapidly as it can around the 
horizon. For example, even in the case of a maximally 
spinning black hole the light-crossing time of the disk in- 
ner edge, taken to be the Innermost Stable Circular Or- 
bit (ISCO), is roughly 1.2 days (accounting for the strong 
gravitational lensing), making large- amplitude intraday 
variability unlikely. 

3. MODELING M87 

There are a number of observational constraints upon 
any model for M87. Foremost among these is the ob- 
served Spectral Energy Distribution (SED), shown in 
Fig. 2, exhibiting a nearly flat, but not inverted, ra- 
dio spectrum and a steep non-thermal decline from mil- 
limeter to infrared wavelengths. It is important to note 
that there is no contemporaneous measurement of M87's 
SED from the radio to the infrared, the points shown in 
Fig. 2 having been collected over nearly two decades. 
The error bars shown represent the instantaneous mea- 
surement errors, and are considerably smaller than the 
observed variability. Further difficulties arise in isolating 
the "nuclear" emission, that arising in the vicinity of the 
black hole, from the emission due to knots in the parscc- 
scale structure of the jet and dust in the nucleus of the 



surrounding galaxy. This latter effect is especially prob- 
lematic between 24 and 100 microns (3 x lO^-lO^GHz), 
where the SED is dominated by nuclear dust emission 
(Perlman et al. 2007). For this reason, we do not at- 
tempt to fit the SED at these wavelengths. 

In addition to the SED, there are a number of 7 mm 
VLB A images, showing the beginnings of the jet on mil- 
harcsecond scales (Junor et al. 1999; Ly et al. 2004; 
Walker et al. 2008). Any viable model for M87 must 
be able to reproduce this extended 7 mm emission, and 
in particular showing the wide jet morphology. 

In Eddington units, M87 is significantly under- 
luminous for an accreting black hole, though consider- 
ably more luminous than a rescaled Sgr A*. This may 
be in part due to the presence of a jet. However, a num- 
ber of authors have found that M87's SED is best fit by a 
two-component model consisting of a Radiatively Ineffi- 
cient Accretion Flow (RIAF)-type disk producing the flat 
radio spectrum, and a jet producing the millimeter, in- 
frared and optical spectrum (Yuan 2000; Reynolds et al. 
1999). The requirement that these models be able to 
qualitatively reproduce the 7 mm VLB A image implies 
that the jet luminosity must be comparable to, if not 
dominate, that from the disk component at this wave- 
length. In the following subsections we describe our mod- 
els for the disk and jet structures, their emission and our 
effort to fit all observational constraints simultaneously. 

3.1. Radiative Transfer 

The primary emission mechanism is synchrotron, aris- 
ing from both thermal and non-thermal electrons. We 
model the emission from the thermal electrons using the 
emissivity described in Yuan et al. (2003), appropriately 
altered to account for relativistic effects (see, e.g., Brod- 
erick & Blandford 2004). Since we perform polarized 
radiative transfer via the entire complement of Stokes pa- 
rameters, we employ the polarization fraction for thermal 
synchrotron as derived in Petrosian & McTiernan (1983). 
In doing so we have implicitly assumed that the emission 
due to thermal electrons is isotropic, but this simplifying 
assumption is not expected to affect our results signifi- 
cantly. 

Regarding the non-thermal electrons, we follow Jones 
& Odell (1977) in assuming a power-law electron distri- 
bution with a cut-off below some minimum Lorentz fac- 
tor. As shown below, the corresponding low-frequency 
cut-off is critical to fitting M87's millimeter spectrum. 
For both thermal and non-thermal electrons, the absorp- 
tion coefficients are determined directly via Kirchoff's 
law (Broderick & Blandford 2004). 

3.2. Disk 

Given the low core luminosity of M87 in Eddington 
units, we employ a radiatively-inefficient accretion flow 
model similar to the one we have previously used for 
Sgr A* (Broderick & Loeb 2006a; Broderick et al. 2008). 
This model is adequate at the low accretion rate of inter- 
est here and makes use of the homologous nature of black 
hole accretion flows. Furthermore, this model incorpo- 
rates, in a limited sense, the presence of outflows, clearly 
appropriate for M87. Unlike Sgr A*, M87's radio spec- 
trum is not inverted, and so the model parameters are 
different. It is important to emphasize that our images 
are largely insensitive to the particular model we adopt 
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for M87's accretion flow, since at millimeter wavelengths 
M87 is dominated by the jet emission. Thus, while we 
include a disk model for completeness, alternative accre- 
tion models produce nearly identical results at millimeter 
wavelengths. 

In our model, the accreting plasma is characterized 
by both thermal and non-thermal populations of elec- 
trons. Both of these are dynamically unimportant, with 
the underlying ions providing the dominant pressure and 
determining the magnetic field strength. The particle 
densities are taken to be broken power laws, with a break 
at some radius which was subsequently adjusted to 
fit the radio spectrum cut-off, making this qualitatively 
similar to the truncated disk solutions discussed by Yuan 
(2000). Outside of R^, the power law indicies of the ther- 
mal particles were taken from Yuan et al. (2003), leaving 
only the normalizations to be adjusted. To fit the radio 
spectral index, we also modified the radial dependence 
of non-thermal particles. 

Specifically, the density of the thermal electrons is 
given by 



nth. = rith.o e 



{r/R^) 
1 



-0.7 



if r > i?b 
otherwise. 



(1) 



where z is the height, R is the cylindrical radius, and 
r = + R^- The thermal electron temperatures is 
given by 

-0.84 



Te=Te,o (r/J?b)" 



(2) 



which is different from the ion temperature, assumed to 
be at the virial value. The corresponding magnetic field 
strength is determined via 



£2 



/? nth- 



6r 



(3) 



and oriented toroidally, where all lengths are measured 
in units of GMj(? . In all models, /? = 10. 

Similarly, the non-thermal electron density is given by 



1 



i{r>Rb 
otherwise, 



(4) 



where the radial power-law index was chosen to roughly 
reproduce the radio spectrum. The associated spectral 
index is 3 (i.e., the spectral flux of the optically thin disk 
scales with frequency v as oc The minimum 

Lorentz factor was set to 100, in approximate agreement 
with Yuan et al. (2003). Finally, the disk velocity was 
assumed to be Kcplerian outside of the ISCO and to 
follow ballistic plunging inside. 

As discussed in §3.4, qualitatively reproducing the 
7mm observations of Sgr A* requires the disk to be- 
come subdominant by this wavelength. We found that 
this requires i?b — 20 GM/c^. The normalizations for 
the thermal electrons are nth,o = 1-23 x IC^cm"^ and 
2^6,0 = 8.1 X 10^ K. The normalization for the non- 
thermal electron density depends upon the particular 
model being considered (see Table 1). For the case 
where the disk is viewed from 40° from the jet-axis, we 
set rinth.o = 6.1 X lO^cm"^, while for all other cases 
JT^nth.o = 3.8 X lO^cm"'^. Note that for all models, the 
non-thermal electron density is substantially lower than 
its thermal counterpart. 



The associated disk spectrum is shown by the red, 
short-dashed line in Fig 2. As required, it becomes sub- 
dominant near 7 mm, and completely negligible by 1 mm. 
For this reason, the particular details of the disk model- 
ing are unlikely to affect the millimeter wavelength im- 
ages of M87's jet. 

3.3. Magnetically Dominated Jet 

While a variety of mechanisms for the formation of 
ultra-relativistic outflows have been proposed, the most 
successful are the magnetically driven models. Despite 
considerable uncertainty regarding their detailed struc- 
ture, magnetic jets share a number of general features, 
which we seek to reproduce here. These include pro- 
viding a qualitatively realistic relationship between the 
acceleration and coUimation, the presence of helical mo- 
tion, the conversion from a poloidally dominated mag- 
netic field to a toroidal magnetic field at large distances, 
and a transition from a jet to a trans-relativistic disk 
wind at large values of the cylindrical radius. Finally, 
any such model must be numerically well behaved every- 
where. 

In practice, all of these conditions could be met by 
appropriating the output of recent jet-formation simula- 
tions (e.g., McKinney & Blandford 2008; Tchekhovskoy 
et al. 2008; Igumens'hchev 2008; McKinney 2006; Haw- 
ley & Krohk 2006; Nishikawa et al. 2005). While this 
approach will undoubtedly prove productive in the fu- 
ture, it is poorly suited to the goals of this investigation. 
Presently, despite being self-consistent, simulations are 
expensive, do not provide sufficiently fine-grained con- 
trols over the jet parameters and would not allow for 
comparisons outside the realm of MHD jets. In ad- 
dition, there is considerable uncertainty regarding the 
mass- loading, critical to the creation of images. For these 
reasons, we opted to produce a qualitatively correct car- 
toon of magnetically dominated jets based upon force- 
free models of Tchekhovskoy et al. (2008), exhibiting 
many of the required properties but providing sufficient 
freedom to explore a variety of jet morphologies that may 
be outside of the purview of existing General Relativistic 
MHD simulations. Here we only summarize the physical 
reasoning and the final structure of the jet, leaving many 
of the details to the Appendix. 

Numerical simulations suggest that outside of the ac- 
cretion disk the plasma is magnetically dominated (see, 
e.g., McKinney & Blandford 2008; McKinney 2006; Haw- 
ley & Krolik 2006). Thus, we model the outfiow as 
stationary, axisymmetric and force-free. In such con- 
figurations, the magnetic field is described by a stream 
function, ip, which must satisfy a 2nd-order partial- 
differential equation. A physically relevant approximate 
solution is 

^ = r2-2«(i_cos6') . (5) 

For < ^ < 1, this solution is within 10% of the 
stream function inferred from existing numerical simu- 
lations (Tchekhovskoy et al. 2008). At large distances, 
R (x z^, and thus ^ defines the collimation rate. For ^ = 
and ^ = 1, the jet is cylindrical and conical, respectively; 
a more typical value is 1/2. 

The structure of the magnetic field and the velocity 
of the outfiow depend upon the angular velocity at the 
field- footprint, 0,{tjj), which is a function solely of tjj. The 



Imaging the Black Hole Silhouette of M87 



5 





-1 -2 -3 

ARA (mas) 

Fig. 3. — Above: 7mm (44 GHz) images for the models listed in 
Table 1 and described in §3.4, convolved with the beam described 
in Walker et al. (2008) (shown in the lower-right). The contours are 
at 0.1, 1, 2, 2.8, 4, 5.7 and 8 mjy/beam (thick lines), and spaced in 
factors of i/2 thereafter (thin lines). In all figures, the jet is oriented 
along the x-axis (East- West), and the black hole is located at the 
origin. The counter-jet is not visible in any of the images because it 
moves away from the observer and lacks the rclativistic beaming of 
the forward jet. These should be compared with the 7 mm image at 
left, produced from the average of 10 images taJsen througout 2007 
and described in Walker et al. (2008) (see also Ly et al. 2004). For 
reference, 1 mas corresponds to roughly 220 Schwarzschild radii. 
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associated four-velocity is 

<=u*j.(l,0,0,f^(V')) 



(6) 



where Up is determined via a normalization condition. 
Near the light cylinder this choice of Up transitions from 
time-like to space-like. There is no physical problem 
with this choice, meaning only that we have chosen a 
pathological observer. Thus, we choose u^. > and set 
a = UpUpfi = ±1. At the light cylinder, this is ill- 
defined. However, as shown below, the plasma structure 
is continuous across this siirfacc, and thus wc simply in- 
terpolate across the light cylinder to obtain its character- 
istics at this position. We define ^{ip) in the equatorial 
plane, choose it to be the Keplerian angular velocity out- 
side of the ISCO, and fix it to the ISCO value inside. This 
boundary is usually what demarcates the "jet" , separat- 
ing it from the surrounding disk wind. 

The stream function also defines the magnetic field 
structure. Within the UF-frame (which we shall denote 
by the subscript F), the poloidal magnetic field 4- vector 
(see the Appendix for the definition in terms of the elec- 
tromagnetic field tensor) is given by 



b"-^ = -Bq 



bl = Bo 



dip 



-BacT- 



sin 61 



-gd9 

^^=So(2-20ar 



cos( 



(7) 

The toroidal magnetic field may be constructed by not- 
ing that hp 0/ti^p must bo a function of i]j alone as well. 
Choosing this function amounts to choosing the current 
distribution within the outflow. However, the condition 
that the poloidal field dominates at large distances gives 



(8) 



where the minus sign arises from the fact that field lines 
are swept back by the motion of the accretion disk. Fi- 
nally, bp is obtained from the identity Uphpfi = 0- 

We assume that the plasma velocity is similar to the 
drift velocity as seen in the frame of a zero-angular mo- 
mentum observer (ZAMO). This is equivalent to assum- 
ing that the plasma is not relativistic at the jet footprint 
within this frame; it gives 



7« 



Pb'^p) 



where 



7 ' 



and P = 



bjrU^r 



(9) 



(10) 



As shown in the Appendix, unlike Up this velocity is al- 
ways time-like. The magnetic field 4-vector in the plasma 
frame is then 

b^=l{b''p-aPhlu'i) . (11) 

All that remains is to define the mass content of the 
outflow. To do this properly would require a full un- 
derstanding of the mass-loading and particle accelera- 
tion in jets, that presently does not exist. Therefore, we 
employ a simple, parametrized prescription based upon 
mass conservation. As shown in the Appendix, the par- 
ticle density n satisfies in Boyer-Lindquist coordinates, 



7n 
u^pb% 



(5" + 5*V^^)=i^(^), 



(12) 



for an arbitrary function of the stream function. Wc 
choose the particular form of F{i})) by defining the jet 
density n on a slice at some height above the disk. Specif- 
ically, we chose the profile to be a Gaussian of the form 

n oc exp ^— r^sin^^/2rfp^, at a height |rcos^| = rfp. 

With u^p and known, this is sufficient to define 
i^(i/)) everywhere. Finally, this was supplemented with a 
cut-off for r < rip, mimicking the loading of the jet above 
the black hole, yielding 



n ■■ 



no 



u'pblFji;) 



(l_e-'-V2.|p) 



(13) 



The shape of the non-thermal electron distribution 
within jet cores is presently the subject of considerable 
debate. Even within the context of power-law models, 
which we employ here, observational estimates of the 
minimum Lorentz factor can vary from 10 to 10^ de- 
pending upon the object, emission model and type of 
observation (see, e.g., Falcke & Biermann 1995; Plan 
et al. 1998; Krawczynski et al. 2002; Jorstad & Marscher 
2004; Sokolov & Marscher 2005; Sauge & Henri 2004; 
Katarzynski et al. 2006; Giebels et al. 2007; Mueller & 
Schwartz 2008; Stccnbrugge et al. 2008). For concrcte- 
ness, similar to the disk wc choose the minimum Lorentz 
factor to be 100, which is broadly consistent with the 
spectral characteristics of M87, as we will see in the 
following section. In contrast to the disk, however, wc 
choose a spectral index between 1.0 and 1.13, depending 
upon the particular model in Table 1. 

In summary, the parameters that control the jet are the 
coUimation rate, ^, the footprint size, rfp, the magnetic 
field and density normalization. Bo and no- The latter 
two are set by fitting M87's spectrum, leaving the first 
two unconstrained. 

3.4. Fitting Existing Constraints 

The two empirical constraints that we considered were 
M87's SED and the 7 mm VLBI images. These are 
not completely independent since, as alluded to earlier, 
the jet at 7 mm must be sufficiently bright to produce 
the extended emission that is observed in the VLBI im- 
ages. Conversely, this requires that the disk become sub- 
dominant, and thus is the motivation behind the trun- 
cated disk. In addition, the cut-off in the spectrum due 
to the minimum Lorentz factor of the non-thermal jet 
electrons is critical to forcing the jet component of the 
spectrum to turn over near 1 mm, ensuring that it does 
not over-produce the flux at longer wavelengths. 

The particular normalizations for the electron densities 
and magnetic fields in the disk and jet are determined by 
fitting the SED shown in Fig. 2. The disk component is 
shown by the red, short-dashed line and falls off rapidly 
near 7 mm. The jet component for our conical jet model 
(MO, Table 1) is shown by the blue, long-dashed line. 
This fits the millimeter-optical spectrum, and is dom- 
inated by the disk spectrum above 7 mm. Within this 
model there; is sonic; ck;tailc;d balancing of the two com- 
ponents, which is presumably associated with the physics 
of jet formation. 

Despite treating the disk and jet components indepen- 
dently, making no attempt to merge the two, we find 
that the jet magnetic field and the magnetic field in the 
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Fig. 4. — Density, velocity and magnetic structure for the jet models listed in Table 1. For each, the density is shown in the left panel 
with thin/thick contours for log(n/nniax) in steps of 0.1/1. The velocity is shown in the center panel, ii' = ^gau^ /u* , with contours of 
the asymptotic Lorentz factor (—ut) in steps of 1 (unit asymptotic Lorentz factor is shown by the lowest contour that extends across the 
entire range in R). The magnetic field structure is shown in the right panel. The contours denote the magnitude of the magnetic field in 
the plasma frame, and are distributed in the same fashion as the density contours. Due to the considerable change in magnitude over the 

range shown, the magnetic field direction is shown by the vector field, where fe' = yJguV" ■ For both the velocity and magnetic field plots, 
the toroidal and poloidal vectors are normalized relative to each other such that they are equal when the toroidal circle has a diameter 
equal to the length of the poloidal vector. In all panels, the magnetic field surface passing through the ISCO is shown by the dashed red 
line. 



TABLE 1 

Jet model parameters 



Model 


a(M) 


e 


rfp(GAf/c2) 




MO 


0.998 


25° 


10 


1/2 


Ml 


0.998 


25° 


20 


1/2 


M2 





25° 


10 


1/2 


M3 


0.998 


40° 


10 


1/2 


M4 


0.998 


25° 


10 


5/8 


M5 


0.998 


25° 


10 


3/8 



inner part of the disk are quite similar, differing by fac- 
tors of a few. In contrast, the jet density is typically 
considerably lower than disk densities, being an order of 
magnitude smaller than that of the non-thermal electrons 
in the disk. However, this disparity is consistent with jet 
simulations which have generally shown that it is quite 
difficult for disk material to diffuse onto the ordered field 
lines within the jet (McKinney & Gammie 2004; Komis- 
sarov 2005; De ViUiers et al. 2005; McKinney 2006). 

We consider six jet-disk models altogether, the param- 
eters of which are listed in Table 1. These explore the 
dependence of the jet structure and images upon the vari- 
ation of a single jet parameter at a time. The spectra for 
all models are nearly indistinguishable from that shown 



in Fig. 2. The 7 mm image, convolved with the beam 
reported in Walker et al. (2008) and Ly et al. (2004) (see 
also Junor et al. 1999) is shown for each model in Fig. 
3. Generally we find an acceptable qualitative agree- 
ment, which is not trivial given the dubious extension 
of the force-free jet model to these angular scales (corre- 
sponding to many hundreds of Schwarzchild radii). We 
note, however, that the brightness profile of our canoni- 
cal jet model decreases much more rapidly (oc z~'^) than 
that observed by Walker et al. (2008) (cx z~^-^), imply- 
ing that corrections to our simplistic model are required 
far from the black hole. Indeed, because we don't model 
the instabilities and subsequent internal shocks noted in 
simulations, which ultimately redistribute the jet's mag- 
netic and particle energy, our jet continues to accelerate, 
reaching unphysically large Lorentz factors by these an- 
gular scales. Nevertheless, we obtain roughly the correct 
morphology for the single jet orientation we consider. 
East- West. 

The density, velocity and magnetic field structure, 
shown in Fig. 4 for the scale that dominates the mil- 
limeter images, share a number of common features. In 
all, the critical surface defined by that poloidal field line 
which passes through the ISCO, shown by the red dashed 
line, is associated with a qualitative change in the out- 
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flow. Following McKinncy (2006), by "jet" we will refer 
the region inside this surface, and by "wind" we mean 
the region outside. It is along this surface that the high- 
est Lorentz factors are reached and the velocity is most 
poloidal (i.e., the velocity pitch angle is largest). As we 
move away from this surface, both within the jet and the 
wind, the toroidal velocity component increases. Con- 
versely, the magnetic field is most toroidal (i.e., the mag- 
netic pitch angle is smallest) along this surface, with an 
increasing poloidal component on both sides. Finally, 
the density generically reaches a maximum along this 
surface. 

This behavior is not unexpected, deriving directly from 
the velocity of the magnetic field footprints. Outside 
of the ISCO, the footprint velocity is monotonically de- 
creasing with radius (since fl is set to the Keplerian 
value). Inside of the ISCO fl is held fixed and thus the 
velocity is increasing with radius. Thus, the critical mag- 
netic field surface corresponds to the maximum footprint 
velocity, with the jet and wind structure on either side 
roughly mirroring each other. 

Next we discuss the individual model features sepa- 
rately. 

3.4.1. MO 

MO represents our canonical jet model, having pa- 
rameters that arc are broadly consistent with numeri- 
cal simulations and direct observations of M87. In this 
case the black hole is rapidly rotating, has a footprint 
of rfp = IQGM/c? and coUimation index of ^ = 1/2. 
Spectra and images are produced by viewing the jet at 
25° from the jet/spin axis, as inferred from observations 
of super-luminal knots (Heinz & Begelman 1997). The 
asymptotic Lorentz factor reaches 5 within the region 
shown in Fig. 4, and is comparable to that inferred in 
M87 (Heinz & Begelman 1997). This model easily pro- 
duces the extended 7 mm emission, with the slight devi- 
ation from purely East- West orientation arising due to 
the helical motion of plasma in the jet. 

3.4.2. Ml 

Ml is identical to MO, save for the increased footprint 
size, rfp = 20. This may be the case, e.g., if black 
hole spin doc;s not play the central (anticipated) role 
in launching the jet. The velocity and magnetic field 
structures are identical to those of model MO, the only 
difference being the width of the density distribution. 
This has consequences for the normalizations of the den- 
sity and magnetic field, with the larger geometric area 
forcing these to be substantially smaller in magnitude. 
Again, this model is qualitatively able to reproduce the 
7 mm observations, perhaps even better than model MO. 

3.4.3. M2 

M2 corresponds to a jet from a non-rotating black hole. 
The primary consequence of setting a = is to increase 
the ISCO radius, thus substantially decreasing the max- 
imum footprint velocity and widening the critical mag- 
netic surface. Thus, even with rip = 10 and ^ = 1/2, 
as in MO, the jet appears considerably broader, and is 
accelerated much more slowly. In the region shown by 
Fig. 4, the maximum asymptotic Lorentz factor reached 
is only 2. Nevertheless, M2 is capable of producing the 
extended 7 mm emission that is observed. 



3.4.4. M3 

The structure of the jet in M3 is identical to MO, the 
only difference being the viewing angle. For M3, this 
is chosen to be 40° from the jet/spin axis. As a conse- 
quence, the density and magnetic normalizations of the 
jet and the underlying disk needed to be adjusted. Oth- 
erwise, there is no difference in the modelling. 

3.4.5. M4 

M4 is a moderately less collimated jet, with ^ = 3/8. 
This relatively small change in the collimation rate has 
dramatic consequences for the density and magnetic field 
structure, both falling off much more rapidly in this case. 
The shape of the asymptotic Lorentz factor contours arc 
somewhat different from those associated with MO, ex- 
hibiting a more gradual transition between the jet and 
wind regions. Nevertheless, the maximum asymptotic 
Lorentz factor reached is 5, similar to MO. In this case, 
the 7 mm image is considerably different than that of the 
previous models, appearing more compact due primarily 
to the rapid density fall-off with radius. 

3.4.6. M5 

Finally, M5 corresponds to a moderately more colli- 
mated jet, with ^ = 5/8. Again this results in a substan- 
tial change in the density and magnetic field structure, 
this time narrowing the jet. Unlike M4, the asymptotic 
Lorentz factor is also significantly altered, reaching only 
3 in the region shown. Thus, the more rapid collimation 
results in a more slowly accelerating, narrow jet. This 
overproduces the extent of the 7 mm emission, and has 
difficulty producing the observed breadth. 

4. JET IMAGES 
4.1. 1.3 mm 

Images of M87 at 1.3 mm (230 GHz) arc presented in 
Fig. 5. As with images of accretion flows (Broderick & 
Loeb 2006a; Broderick et al. 2008), the outflow images 
are strongly influenced by relativistic beaming, boosting 
and strong gravitational lensing. In all images the jet 
is lopsided, with the emission dominated by the side in 
which the helically moving outflow is approaching the ob- 
server. However, whereas for accretion flows this results 
in a bright crescent, for our outflow model this produces 
a bifurcated image. This is clearly seen in the 1.3 mm im- 
age of MO, showing bright diagonal and vertical branches, 
which can be understood in terms of the generic velocity 
structure of the jet. As discussed in §3.4, the velocity 
is both largest and nearly radial along the critical mag- 
netic field surface that divides the jet from the wind. 
Here, relativistic beaming dims the emission outside of 
a small angle, producing a flux deflcit. Both inside and 
outside this surface, the toroidal component of the veloc- 
ity increases, and the velocity is aligned along the line of 
sight, producing bright regions. The bright vertical and 
diagonal branches are due to emission from the jet and 
wind, respectively. 

Even at a wavelength of 1.3 mm, the jet is typically not 
completely optically thin. This is not unexpected given 
that the jet spectrum is just turning over at 1.3 mm. This 
is most pronounced near the jet base, which obscures a 
portion of the black hole horizon. Nevertheless, in all im- 
ages the silhouette cast by the black hole upon the lensed 
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Fig. 5. — 1.3 mm (230 GHz) images for the models listed in Table 1 and described in §3.4. The brightness scale is logarithmic, showing a 
dynamic range of 128, similar to the existing 7 mm images. In all images the projected jet axis is parallel to the vertical axis. The ellipses 
in the lower left provide an estimate of the beam size for arrays consisting of the North American telescopes only (blue), including the 
European telescopes & the LMT (green), and all telescopes (yellow). 



image of the base of the counter-jet''^ is clearly visible. 
This suggests that imaging a black hole silhouette may 
be more easily done in M87 than in Sgr A*. However, at 
1.3 mm, the optimal beam, obtained by phasing together 
an array including all existing North American and Eu- 
ropean telescopes, the LMT and ALMA, is comparable 
in extent to the silhouette. 

For the purpose of informing jet-formation modelling, 
it is fortuitous that each of the models is clearly distin- 
guishable. The wider footprint of Ml results in a corre- 
spondingly broader brightness distribution. The lower 
velocities associated with lack of spin (M2) produces 
a featureless, asymmetric image with a single branch 
(though completely optically thin). Changing the incli- 
nation (M3) produces a broader base and narrower jet 
branch, though this case is most degenerate with the im- 
age from the canonical model (MO). It is not obvious 
whether images of the quiescent jet emission alone will 
distinguish these two cases. On the other hand, alter- 
ing the coUimation index even a small amount produces 
drastically different millimeter images (M4 & M5). Thus, 

* In this case, we are seeing the counter-jet near to the black hole, 
where it has not yet accelerated appreciably and thus relativistic 
beaming is not yet effective. Far from the black hole the jet is 
rapidly moving away from the observer and beaming renders it 
invisible, as seen in Fig 3. Should the jet accelerate much more 
rapidly than modeled here, or the counter-jet not be present, the 
silhouette will be much weaker. 



short-wavelength imaging promises to severely constrain 
the jet footprint size, the black hole spin and the coUi- 
mation index, all providing critical observational input 
into existing efforts to model the structure of MHD jets. 

4.2. 0.87 mm 

All of the features seen at 1.3 mm are also present at 
a wavelength of 0.87mm (345 GHz). The two sets of im- 
ages are subtly different primarily due to the decreased 
optical depth. Despite this, there are still optically thick 
regions for all models except M4. Again it is possible 
to clearly distinguish footprint size, black hole spin and 
collimation model. At 0.87 mm, the jet orientation is 
more easily discerned, though it remains unclear how ac- 
curately it may be determined. 

Instead, most important is the reduction in the beam 
size relative to the horizon scale. At submillimeter wave- 
lengths the black hole silhouette is now larger than the 
beam, implying that it is possible to clearly identify it in 
these images. However, comparison of the M4 image and 
the others suggests that it will not be possible to mea- 
sure black hole spin with the silhouette detection alone, 
whose shape and size is similar for both high and low 
spins. 

This is not unexpected since the primary effect of spin 
is to shift the silhouette. Unlike the horizon, the appar- 
ent silhouette size is insensitve to the spin magnitude. 
This can be easily see by considering critical orbits in 
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Fig. 6. — 0.87 mm (345 GHz) images for the models listed in Table 1 and described in §3.4. The brightness scale is logarithmic, showing 
a dynamic range of 128, similar to the existing 7 mm images. In all images the projected jet axis is parallel to the vertical axis. The ellipses 
in the lower left provide an estimate of the beam size for arrays consisting of the North American telescopes only (blue), including the 
European telescopes & the LMT (green), and all telescopes (yellow). 



the equatorial plane. The black hole silhouette is pro- 
duced by geodesies which terminate on the horizon. For 
a non-rotating black hole, these correspond to impact pa- 
rameters less than \/27GM/c^ ~ 5.2GM/c^, and thus a 
silhouette diameter of roughly 10.4 GM/c^. For a rapidly 
rotating black hole, the critical impact parameter for pro- 
grade geodesies decreases. However, this is compensated 
for by the increase in the critical impact parameter of 
retrograde geodesies. In the case of a maximally rotat- 
ing black hole (a = 1) these are 2 GM/c^ and 7GM/c^ 
respectively. Thus the associated silhouette diameter has 
shrunk only to 9GM/c^, though it is shifted by roughly 
2.2 GAf/c^ (see, e.g., Chandrasekhar 1992). As a result, 
efforts to discern the black hole spin by careful measure- 
ments of the silhouette size and shape alone are likely to 
be compromised by the rather large astrophysical uncer- 
tainties associated with the mass and distance to M87's 
black hole, as well as the resolution of future experiments. 

5. JET POLARIZATION 

Polarization maps for the six models we consider are 
shown in Fig. 7. The polarized flux can be compara- 
ble to the total flux, with the polarization fraction being 
largest where the flux is greatest. As with earlier radio 
polarization studies, the polarization is indicative of the 
magnetic field structure, and in particular, the relation- 
ship between the plasma density and the magnetic field. 



The maximum polarization fraction in all images is 
roughly 75%, is strongly correlated with the intensity, 
and thus should be relatively easily detected. The domi- 
nant magnetic field orientation is easily determined from 
the polarization maps in the usual way: assuming that 
the projected magnetic field direction is orthogonal to the 
polarization vectors. Unlike in some quasar jets (Lyu- 
tikov et al. 2005), the polarization angles here are not 
significantly altered by the relativistic motion. This is 
because of: (i) the moderate Lorentz factors in the re- 
gions contributing to the millimeter emission, and (ii) 
the orientation of the jet near the observer's line of sight. 
Most importantly, polarization measurements appear to 
provide a method by which the structure of the magnetic 
field lines may be compared directly with the coUimation 
of the emitting plasma. 

Since in our jet models the emitting electrons are rel- 
ativistic, it is not possible to produce a robust estimate 
of the rotation measure (RM) from thermal electrons, 
and thus infer the importance of Faraday rotation for 
polarization measurements in M87. If there is a thermal 
population of electrons with a density comparable to the 
non-thermal component, we estimate a RM of roughly 
lO^rad m~^, which produces a rotation in the plane of 
polarization at 1 mm of 60°. In this case, Faraday depo- 
larization is unlikely to be a significant problem, and it 
will be possible to reconstruct the unrotated polarization 
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Fig. 7. — 0.87mm (345 GHz) polarization maps of the models in Table 1, superimposed upon logarithmic contours of the intensity (factors 
of 2). The size of the tick marks are indicative of the polarization fraction, with the maximum polarization typically being 75%. The 
polarization vectors have not been Faraday rotated. The ellipses in the lower left provide an estimate of the beam size for arrays consisting 
of the North American telescopes only (blue), including the European telescopes & the LMT (green), and all telescopes (yellow). 



map. 

This should be tempered, however, by measurements 
of the rotation measure within M87's jet at longer wave- 
lengths. At 20 mas from M87's core, the RM has been 
found to range from — 4 x lO'^rad m~^ to lO'^rad m~^, 
depending upon the location observed (Zavala & Tay- 
lor 2002). If these are not enhanced locally, due per- 
haps to a substantial over density of cold electrons within 
the knots being observed, they imply a rotation measure 
orders of magnitude higher near the black hole. That 
is, if n cx and B cx r~^, as expected at large dis- 
tances, RM cx nBr oc r~^. Extrapolating from 20 mas 
to 200 /ias implies a nuclear RM of 10* rad m~^, two or- 
ders of magnitude larger than our estimate. In this case, 
band-depolarization will still not be important for the 
largest bandwidths likely to be accessible (a few GHz). 
However, in situ Faraday depolarization will play a role, 
possibly attenuating the polarized flux considerably, de- 
pending upon geometric considerations. As such, the 
absence of polarization itself will constitute an interest- 
ing and potentially important constraint upon the non- 
relativistic electron population within M87's outflow. 

6. CONCLUSIONS 

M87 provides a promising second target for the emerg- 
ing millimeter and submillimeter VLBI capability. Its 
presence in the Northern sky simplifies its observation 
and results in better baseline coverage than available for 
Sgr A*. In addition, its large black- hole mass, and corre- 



spondingly long dynamical timescale, makes possible the 
use of Earth aperture synthesis, even during periods of 
substantial variability. 

Despite being dominated by a jet at millimeter wave- 
lengths, the black hole silhouette, cast against the base 
of the counter-jet, is still clearly present, and may just be 
resolved at 0.87mm (345 GHz). This requires that both 
the jet and the counter-jet are moving slowly near the 
black hole, and thus relativistic beaming does not pre- 
clude its observation. If, however, the counter-jet is ei- 
ther relativistic even at its base or simply not present the 
silhouette will not be detectable. Nonetheless, obtaining 
a sub-mm image of M87 appears to be more promising 
than efforts to image the silhouette of Sgr A* due to the 
lower opacity of the surrounding synchrotron emitting 
gas. Unfortunately, it is unlikely that the silhouette alone 
will be sufficient to constrain the spin without additional 
information. Better constraints may originate from ex- 
ploring the effect of spin upon the jet image (via a change 
in the angular velocity profile of the jet footprint) or by 
carefully studying variability of hot spots in the accre- 
tion flow and jetted outflow. The latter possibility is fea- 
sible owing to the low optical depth to synchrotron self- 
absorption and the long dynamical timescales in M87. 

Direct imaging of M87's jet on horizon scales will neces- 
sarily provide a wealth of information crucial to modern 
numerical efforts to model jet formation and propaga- 
tion. In particular, such observations will be sensitive 
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to the jet footprint size, coUimation rate, and (if the jet Ultimately, such measurements will associate the black 

is truly a magnetically dominated outflow) black hole hole demographics measured by LISA with the existing, 

spin. Polarized imaging will further shed light upon the extensively studied set of known quasars. 

structure of the jet's magnetic field and its relation to 

the emitting gas. Less well constrained will be the jet 

orientation. 

Finally, comparisons of M87 and Sgr A* will provide Wc would like to thank Bill Junor, Craig Walker and 

the first opportunity to relate the gravitational prop- Chun Ly for graciously allowing us to use their VLBA 

erties of quasars, about which much is known, to the images of M87. This work was supported in part by 

lower-mass, dormant black holes that LISA will access. Harvard University funds. 

APPENDIX 

FORCE-FREE JET IN COVARIANT NOTATION 

In order to perform the radiative transfer through the region around the M87 black hole we must characterize the 
surrounding plasma. This is done in terms of a proper density, n, a four- velocity, u'^, and the magnetic field within the 
plasma frame b^^ = Ui,*F^'^, where *F^'^ is the dual of the standard electromagnetic field tensor. Specifying that the 
magnetic field be force-free simplifies the definition of these quantities substantially. However, usually this is presented 
in a 3-dimcnsional vectors instead of the fully covariant form we require. In this appendix we review some relevant 
theorems and determine the jet structure in both notations. 

Vector Force-Free Theorems 

The "Force Free" condition is 

pE-hjxB = 0, (Al) 

i.e., the Lorentz force vanishes. A direct consequence is that E and B are orthogonal: 

B • (pE-hj X B) = pB • E = E-B = 0. (A2) 

This implies that there exists a frame in which one of these fields vanishes. For a force-free solution to exist we require 
that in this frame E vanishes, otherwise there would be no way to satisfy (Al) generally. Thus, the Lorentz scalar 
^2 _ ^2 jjjyg^ non-negative. The velocity of this frame, called the "drift velocity" is not unique. However, we can 
define a class of such velocities: 

E X B 

VF = — + ^ E = VF X B . (A3) 

Choosing /? determines the frame in which we seek to define the drift velocity. 

Maxwell's equations may be combined with the force-free condition to produce equations in which only the electro- 
magnetic fields appear. For stationary solutions these give 

V-B = 0, VxE = and Bx(VxB)=EV-E, (A4) 

where the last of these arises directly from the force-free condition itself. This may be rewritten in terms of an unknown 
scalar function, a, as 

V X B = aB-Fv^V -E, (A5) 
where a and /3V • E are degenerate. We shall refer to this second equation as the Force-Free equation. 

Conserved Quantities 

In axisymmetry, V x E = implies that there can be no toroidal component to E, and thus we may choose vf to 
be purely toroidal as well. There arc a number of immediate consequences. Firstly, 

B • Va = V- aB = V • (V X B- vfV • E) = 0, (A6) 

and a is constant along poloidal magnetic field surfaces. Secondly, inspecting the poloidal components of the Force-Free 
equation (in spherical coordinates for convenience): 

\/hB'' = -deBA, and \/hB^ = --drBs , (A7) 
a a 

where h is the determinant of the 3-metric, thus is constant along poloidal field surfaces as well. Coupled with the 
previous result for a, this implies that there exists a "stream function", tp = J a~^dB^, which satisfies: 

VhB'' = d0^ and VhB^ = -drip . (A8) 

Then, we can define Bff,{ip) by its value on the equatorial plane (or any other convenient surface crossing every poloidal 

field line), which implicitly defines a{tjj) = dip/dB^. Thirdly, defining fl = Vp, this implies: 

E = f70xB ^ VxE = 0B-Vf7 = O, (A9) 

and thus Q is constant along poloidal field surfaces as well. Note that for this choice of vf, there can be light- 
cylinders, a region where vp > c. This does not mean that the Force Free condition has broken down. Rather, as 
we shall explicitly see in the jet models discussed below, the choice of vf may simply be inappropriate. Nevertheless, 
the relationships between E and B are unchanged, and we may construct solutions that are continuous in the relevant 
observable quantities. 
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Determining the Stream Function 
The stream function itself is determined from the phi-component of the Force- Free equation: 



+ 9rr {deV4,) 



If we recast this relation in cylindrical coordinates, {R, z^cj)), it reduces to 



(1 - R^n^) dl^ + (1 - R^Q^) dli, - (1 + R^^f) -Or^ - flR^— (dR^y + + B^-^ = , (All) 



+ ^B'^^ = 0. (AlO) 
d^ 



dtp 



which, after trivial rcnamings, is identical to Eq. (13) in Narayan et al. (2007) (and also Okamoto 1974). All that 
remains is to specify, B^{ip) and solve for ip. 

Vector Jet Structure 

Here wo use the various insights gained in the previous section to write down a qualitative jet model in covariant 
form. It is important to note that in both cases we are attempting only a qualitatively correct "solution" , paying 
specific attention only to ensuring that the magnetic field geometry and the jet velocity do not become unphysical 
anywhere. 

Following Tchckhovskoy ct al. (2008), we begin by positing a class of stream functions based upon solutions to the 
= stream-function equation: 

V' = rP(l-cos6') . (A12) 
From this, we may immediately construct the poloidal magnetic field structure: 

B^ = rP-^ and = -pr^-^ ^ ^ ^ = -pr^"^ tan - . (A13) 

smw 2 

Since the toroidal component of the magnetic field scales as r~^, far from the jet source it becomes electromagnetically 
dominated. Thus, B'^B^ ~ £^ = R^Q^B^Br, and B^ = —R^flB^, where the negative sign arises from the fact that 
the field lines are swept back (recall that R = rsinO is the cylindrical radius). That is, 

B^ ~ -r^ sin^ fl r^'^ ~ -20^ • (A14) 

Therefore, since Btf, is conserved along magnetic field lines, 

B^ = ^ = -M^i = -20(^).-^^^ . (A15) 
i?2 r^s\n^e sine ^ ^ 

Thus, given ^l{tp) (which we shall determine below), we now have the magnetic field value B everywhere. 

Plasma Velocity 

Because the plasma is non-relativistic initially, and the Bernoulli constant is conserved along the plasma flow, the 
plasma velocity is very similar to the drift velocity in the lab frame, and thus it is suflicient to compute the latter. 
The electric field in the lab frame is simply E = il,cj) x B. Therefore, 

v=^ = -^0+^B (A16) 
Plasma Density 

Finally, we must determine the electron/pair density in the jet. We do this assuming that the continuity equation 
holds, i.e., V • nv = 0. This is a dubious assumption for two reasons. First, if the mass-loading of the jet occurs over 
an extended region (e.g., if it is due to photon annihilation), then continuity fails locally. Second, and more severe, 
is that we are interested only in the non-thermal tail of the electron population. Thus, we might expect the electron 
acceleration process and cooling to leave an imprint upon the non-thermal electron distribution (which is completely 
ignored in our prescription). Nevertheless, in the absence of a more detailed understanding of jet formation, we proceed 
with this fiawed working assumption, which gives, 

V-nfi<^ + V-n^B = B. Vn^ = 0, (A17) 

or nQB^/B'^ is constant along field lines. Since Cl and B^ are already known to be constant along field lines, this gives 

n cx 
finding 

n = B'^F{ip) , (A18) 
everywhere else, for some function F{ijj) that depends upon the jet loading at the footprint. 



n oc B^ along the field line. Thus, we must only specify n/B'^, i.e., the squared Alfven velocity, at the jet footprint. 
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Covariant Force-Free Theorems 

Each of the results in the previous sections can be presented in a fully covariant notation as well. In terms of the 
electromagnetic field tensor, the "Force Free" condition is 

= . (A19) 

This also directly implies an orthogonality relation: 

*F,^.F"ri. = -^*F^.F'^''j, = ^ *F^-F^, = 0, (A20) 

where we used the general relation *A^'' A^v = *A°"' Aai,5^ / 4 for all antisymmetric 2nd-rank tensors, A^'^ (which may 
be derived directly from the definition of the dual). This is not particularly surprising given that the Lorentz invariant 

*pnvp^^ = 2E • B, where E and B are determined in the locally flat coordinate patch. 

Again, wo require that it is possible to transform into a frame in which the electric field vanishes, which corresponds 
to the condition that -F^'^-F^y > 0. As before, the velocity of this frame is not miiquely defined, nor must it be physical. 
However, if one such velocity is u^, then we require that the 4- vector that is coincident with E in this frame vanishes, 
i.e., 

e'^ oc u^Ft = . (A21) 
Thus, the acceptable zt^ comprise the null-space of F'^ . Generally, given any u^^ we can write F^^" and *F^'' as 

F"'^ = ?i^e'' - zi-^e^ + e^'''"'\iah8 and *F^'^ = u^'b'' - u^bf" + e^'^'^'^Uc^ep , (A22) 

where e^'^"'^ is the Levi-Civita pseudotensor. If we set = Up, and thus = 0, this implies that we can write these 
as 

p^u ^ £i^'^°'I^UFabF0 and *F'"' = u'^pb^'p - u^'pb'^p . (A23) 

Since we require that there exists a physical frame in which e'^ vanishes, there must be time-like choices for Up. This 
implies that the null space of F™'^ is two-dimensional with one time-like dimension and one space-like dimension (since 
<^ - 0). 

Since we know that *F'^F'^ = from the aforementioned orthogonality condition, the 4-vectors v'^^.^ = *F'^ are 
necessarily components of F'''"s null space. In fact, we can show that this space is completely spanned by the v'^^y 
which are themselves degenerate. This may be proved by inspection, constructing a space-like and time-like basis 
vectors: 6^*Ft = and w^*F't = u^b'^^. 

Regardless of the nature of Up, we may always construct time-like drift velocities. For any time- like vector, r/'*, 

, "-^""-^'^'r ^^^(^u'i+(3b>i) where (3=^MeI1^ and 7= (A24) 

Since the quantity v?p appears in a number of places, let us define a = u\ = ±1. That 77^ be time-like itself is 
sufficient to guarantee that will be time-like as well. To see this, let us decompose ry'^ into a tetrad basis with 



'(0) ~ ""-F' ^(1) ''(2,3) 

\^rj /-|\2q //""\\2r) /i\2r) /Q\2r) /q\2 



b^ and e^^ orthogonal to the rest (and thus necessarily space-like). Then, 



oc ^^el, + r,^')%l, < ry(") e^„) + ^i'^e},) + v^'^ef,) + v^'^ef,) = < . (A25) 

Finally, we have chosen the arbitrary sign of such that it corresponds to forward-propagating particles. To see this 
generally, we again appeal to the tetrad expansion of rj: 

^(.^°^^?o)+.'^'^^,). (A26) 



UpTju 



Since the term in the parentheses is negative and we choose u^r/^ to be negative by convention, this implies that 7 
must have the same sign as a. 

Note that if cx u^, we find = u^. Thus, choosing rj'^ is equivalent to choosing the frame in which we wish to 
define the drift velocity. That these are not all equivalent is a consequence of the ambiguity under boosts along the 
magnetic field. Nevertheless, if we are given a preferred frame (e.g., the lab frame) there is a unique definition, as in 
the 3-D vector computation. By inspection in the locally-flat coordinate frame we can show that this definition of 
is identical to the 3-D vector case. 

We now turn to the various representations of the Force- Free equation. Recalling that Maxwell's equations give us 
Vi/F^'^ = 47rj'', the force-free condition may be trivially written in terms of F^^: 

FtV<.F^'" = 0, or W ^F^" = ab% + Cu% , (A27) 

i.e., once again we find that the vector we desire (in this case j**) is in the null-space of F'^. Now, there are two 
undefined functions, a which is analogous to the 3-D vector case, and C,. As it turns out, we will be able to completely 
determine the magnetic field without ever having to specify Q, which is not surprising given that there are only three 
independent components to b'^. 

^ Note that should we choose a form for u^. that becomes space-hke in some region, 6^, necessarily becomes time-like. 
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Conserved Quantities 

First, let us derive a number of useful relations in axisymmetry. As before, we shall assume that it is possible to put 
the field- line velocity (drift velocity) into a purely toroidal form by an appropriate boost along bp-. 

u^p = u%t^ + u%(j)^ , (A28) 

where and (j)^^ arc the time-like and azimuthal Killing vectors present in all stationary, axisymmetric spacetimes. 
The first thing we note is: 

= bp^W,*F^^ = bp^V. {u^pb-'p - b^pu-'p) = b%b),W,UF^ , (A29) 

and thus the magnetic field must lie upon non-shearing velocity surfaces. This may be simplified by Killing's equation 
and the condition that bp^Up = ^ bpf = —bp^fl where fl = u'^/up-. 

b^pb^pV^up^ = b'^pb-'p (t^d.u'p + </)^a,M^) = bp^b-'p {d^Qu'p - Qd^u'p) = bp^u'pb-'pdM . (A30) 

That is, is constant along magnetic field lines. 
Next we note 

= up^V.*F'"' = up^V. {u'^pb-'p - Wpv^p) = b^pU^pV.up^ + aVub^p . (A31) 
The first term may again be simplified using Killing's equation: 



b'^u'pV^up^ = b'^UpV^ [Uptf, + u'},(f)^^ 

= bp (t^updi,up + (j)^Upd^Up + UpUpVutu + UpUpV^cj)^^ 
= -b^p (u'pU^pVi.t, + u^u^pV 

= b'P (ti,Upd^Up + (j)^Upd^u''p - Up'^^upJ^ (A32) 

= 6^ (uptd^Up + UF^df^v!\ 



'■F 

.t 



(A34) 



= bp {uptd^up + up^d^Qup) 
= b'^F [{upt + up^Q) dfiU^p + up^u\rdij,Q\ 
= ab'^d^ Itiu% . 

Therefore, V^bp = —6^9^ In u^.. However, from the Force-Free equation, we have 

= V^V.F'^'- = V^a6^ + Vi^Cu^p = b%d^a + aV^b^p = b^d^a - al/'pd^ Inu*. = u'pb'^pd^— , (A33) 

Up 

and thus a/u*p is conserved along field lines. 
This now brings us to the stream function. Again, looking at the poloidal parts of the Force-Free equation we find: 

ab"" = s'^''"^d^upabpg = -^=de (upJjpt - uptbpj,) = — ( -7^ ] 

\Up J 

ab^ = s^^^^d^upabpfj = --^=dr [up^bpt - up^bp^) = "^^^^ (^) ' 
where we used 

uptbp^ - up^bpt = uptbp^ + up^bp^^ = — ^ (uptUp + up^uyi = a—^ . (A35) 
Since a/u* is conserved along field lines, we may define a stream function by 

which is very similar to the 3-D vector case, except for the presence of the u^. In terms of this, we have 

:/zli,r -(jg^^ and ^^bi = adrip . (A37) 



Now the problem is again reduced to determining ^. 
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Determining the Stream Function 

At first glance, it may appear that we must also know ( to obtain an equation for tjj. However, we can use the fact 
that bpUp^ = to avoid this requirement: 



{bFt'^F^, - bF^UFt) (drbFe - debFr^ + bFodr 



bpa 



bFr^e 



bpd 



bpe (bptdrUF^ - bi.\i,drUFt) + bpr {bFtdeUF4, - bpadeUFi) 



zdrip + de- 



= ay/^albprbp + bFebp) + <^-i^ 

V /Up 

- bpo (u%drUF4, + UpdrUFtj + bpr [u'^pdeUF^^ + UpdeUFtj 

where the following relations were used, 

bFtdiUFcj, - bF4,diUFt = -0.bF4,diUF<j, - bF(j,diUFt = (u'^diUF4, + u]rdiUFt^ 

w 771 ^ / 



(A38) 



bpftdr 



Note also that 



b% = bFtb% + bp^b^ + bp^¥p + bpeb^p = bp^b^p + bpgb% + ^ - ut,b%) 



(A39) 



(A40) 



and thus, 

dr—^^^j=dri) + dg—j^^=dg'lp - abpg iutdrUpj, + U^pdrUpA 

+ abpr (u^dgup^ + Updgupf 



;i%Mf^^.„. (A41) 

dip \ upt 



This is quite similar to Eq. (AlO), but significant differences remain. Most notably arc the presence of various flow- 
velocity components (the u^'s). When = this simplifies somewhat, though not as completely as in the 3-D vector 
case. The reason is that there may still be non- negligible terms due to gravity (i.e., non- vanishing derivatives of gu)- 
In flat space, however, this gives precisely the same expression as Eq. (AlO). 



Covariant Jet Structure 
Again, let us begin with the the class of stream functions, 

^ = rP(l-cos^) . 

Furthermore, at each location we have u* in terms of x'' and fl{ip): 



Up = \gtt + "^gtrti^ + g^(ji^ 



21-1/2 



(A42) 



(A43) 



where some care has been taken to deal with the possibility of it^ being space-like. We then determine the poloidal 
magnetic field seen by the orbiting observer as 



sin 9 



-ap- 



rP-^l - cose) 



(A44) 



p , f^Pt 



bpf = -bp^Q. 



All of these functions are well behaved in the region of interest. Note, however, that until we settle upon the plasma 
velocity, these results do not provide the magnetic field in the jet frame. 

Plasma Velocity 

Again, we assume that the plasma velocity is similar to the drift velocity as measured in some lab frame. Based 
upon the physical notion that zero-angular momentum gas should not see an acceleration along the jet (and guided 
by the requirement that we have a meaningful time- like vector), we choose 



1 



F, 0,0,0 



(A45) 
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As a result, 



7 {up + (3b p) where (3 = 



hi 



p^p 



and 7 is given by Eq. (A24). In this frame, the magnetic field is 

6" = u^*F'''' = -u^a {u'^ph'^p - b^vTp) = 7 (6^ - ahl(3u^p) , 
where the additional factor of — tr is necessary to keep *F^" continuous across the light cylinder^. 



(A46) 



(A47) 



Plasma Density 

Again, we set the density using the continuity equation. This gives 

V^nui" = V^n7 « + ^h%) = jpriV^b^p + b^pd^jpn = u'p\fpd^ 



7/Jn 



= 0. 



(A48) 



Thus, similar to the 3-D vector case, some combination of the magnetic field, VL and the density is conserved along 
field lines. In Boyer-Lindquist coordinates, we can simplify (3 somewhat to bring this expression more in line with the 
3-D vector case: 



The first two non- trivial factors, bp^/u^p and O, are conserved along field lines, which is very similar to the 3-D vector 
case. Therefore, 



b%u*p b%u*pUFt 



b'^pd^ 



= 0, 



is conserved along field lines. That is, the density is given by 



^p^p 



(A50) 



(A51) 



given the arbitrary function -F(V')- This is chosen as described in the main text. 

Parameters & Summary 

As before we must set three parameters (the power-law index p and the density & magnetic field normalizations) 
and two functions of ip {i^{tp) and F{tp)). We describe how these are chosen in the main text. Given these, in practice 
the procedure for computing the magnetic field, plasma velocity and plasma density at a point x^^ is: 

1. Settle upon p, fi(V') and -F(V')- 

2. Construct bpix^) and u^p{x'^). 

3. Prom these, construct b'^{x^), u'^{x^) and n{x'^). 
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